import numpy as np
import matplotlib.pyplot as plt
import cv2
#初始化显示窗体：激活状态分布MapM和能量分布MapQ
cv2.namedWindow('LifeWorld_MapM', cv2.WINDOW_NORMAL)
cv2.namedWindow('LifeWorld_MapQ', cv2.WINDOW_NORMAL)
#初始化显示信息熵和生命结构数的坐标系
plt.ion() 
figure,ax=plt.subplots(2,1)
plt.suptitle("Entropy And CountOfLife")

rg = np.random.default_rng()

#初始化激活和能量分布
def initMap(w,h,r):
    mapm=np.zeros((w, h), dtype=np.uint8)
    mapq=np.zeros((w, h))
    #随机初始化激活状态
    mapm[rg.random((w,h))>r]=1
    #初始化能量分布
    sumq=np.log(3)/2*w*h
    mapq[0,0]=sumq
    return mapm,mapq

#初始化设定生命结构
def initLife(w,h):
    #0/1：生命结构中必须为该值，2：生物可控制位置，3：任意值
    struct=[[3,2,2,2,3],
            [2,0,1,0,2],
            [2,1,0,1,2],
            [2,0,1,0,2],
            [3,2,2,2,3]]
    life=[[],[],[],[]]
    for i,row in enumerate(struct):
        for j,val in enumerate(row):
            life[val].append((i,j))
    mapctrl=np.zeros((w,h))
    return life, mapctrl

#通用工具：循环平移
def shift(mapm,i,j):
    mapm=np.vstack((mapm[i:,:],mapm[0:i,:]))
    mapm=np.hstack((mapm[:,j:],mapm[:,0:j]))
    return mapm

#计算下个状态的激活状态和能量分布(激活状态，能量分布，热导率)
def nextState(mapm, mapq, k):
    mapmstack=[]
    maptstack=[]
    mapt=np.zeros_like(mapq)
    mapt[mapm==1]=mapq[mapm==1]
    for i,j in ((-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)):
        mapmstack.append(shift(mapm,i,j))
        maptstack.append(shift(mapt,i,j))
    mapmstack,maptstack=np.dstack(mapmstack),np.dstack(maptstack)
    tm=np.sum(mapmstack,axis=2)                         #激活邻域数
    #熵增规则，跟新能量分布
    tq=np.sum(maptstack,axis=2)                         #激活邻域能量和
    mapq[mapm==1]+=(tq[mapm==1]-(mapq*tm)[mapm==1])*k   #激活区域间热传导
    #生命游戏规则，更新激活态分布
    mapm[tm!=2]=0
    mapm[tm==3]=1
    #概率规则，翻转激活状态
    mapp=(2/(np.exp(-2*mapq)+1)-1)                      #翻转状态概率分布
    mapt=rg.random(mapm.shape)<mapp
    mapm[mapt]=np.uint8(np.logical_not(mapm[mapt]))
    return mapm,mapq,mapp

#按比例显示(激活状态，能量分布，可控制区域，缩放比例)
def showMap(mapm, mapp, mapctrl, c):
    R=mapm*255
    G,B=R.copy(),R.copy()
    G[np.logical_and(mapctrl,mapm==1)]=0    #粉色:可控细胞当前激活状态1
    G[np.logical_and(mapctrl,mapm==0)]=255  #绿色:可控细胞当前激活状态0
    imgm=cv2.resize(np.dstack((B,G,R)),None,fx=c,fy=c,interpolation=cv2.INTER_NEAREST)

    imgq=np.uint8(mapp*255)
    imgq=cv2.resize(imgq,None,fx=c,fy=c,interpolation=cv2.INTER_NEAREST)
    
    cv2.imshow('LifeWorld_MapM',imgm)
    cv2.imshow('LifeWorld_MapQ',imgq)

#计算信息熵
def entropy(mapp):
    w,h=mapp.shape
    Emax=2*w*h*0.5*np.log(0.5)
    mapp=mapp[mapp>0]
    mapp=mapp[mapp<1]
    E=(np.sum(mapp*np.log(mapp))+np.sum((1-mapp)*np.log(1-mapp)))/Emax
    return E

#识别生命结构和可控制区域(激活状态，能量分布，生命结构)
def identify(mapm, mapq, life):
    mapt=mapm.copy()
    mapt[mapq==0]=0
    zeros=np.zeros_like(mapt)
    ones=np.ones_like(mapt)
    for i,j in life[0]:
        zeros = np.logical_or(zeros,shift(mapm,i,j))
    for i,j in life[1]:
        ones  = np.logical_and(ones,shift(mapt,i,j))
    maplife=np.logical_and(np.logical_not(zeros),ones)

    mapctrl=np.zeros_like(mapt)
    for i,j in life[2]:
        mapctrl=np.logical_or(mapctrl,shift(maplife, -i, -j))
    return maplife, mapctrl

#运行模拟(宽度，高度，热导率，初始激活比例，显示比，间歇时间，最大迭代数)
def runWorld(w, h, k=.1, r=.86, c=4, wait=0, maxi=10000):
    mapm, mapq   =initMap(w,h,r)
    life, mapctrl=initLife(w,h)
    #仅可控制区域可点击，控制反转状态
    def mouse_callback(event, x, y, flags, param):
        if event==cv2.EVENT_LBUTTONDOWN:
            y,x=x//c,y//c
            if x<w and y<h and mapctrl[x,y]:
                mapm[x,y]=0 if mapm[x,y]==1 else 1
                showMap(mapm,mapp,mapctrl,c)
    cv2.setMouseCallback('LifeWorld_MapM', mouse_callback)
    
    for i in range(maxi):
        mapm,mapq,mapp  =nextState(mapm, mapq, k)
        maplife,mapctrl =identify(mapm,mapq,life)
        #显示信息熵和生命结构数
        E=entropy(mapp)
        C=np.sum(maplife)
        ax[0].plot(i,E,'.')
        ax[1].plot(i,C,'.')
        plt.pause(0.01)
        #显示激活状态和能量分布
        showMap(mapm,mapp,mapctrl,c)
        #暂停的情况
        if E>0.9:
            wait=0
        elif C>0:
            wait=0
        key=cv2.waitKey(wait)&0xFF
        if key==ord('q'):   #退出
            break
        elif key==ord('p'): #暂停
            wait=1 if wait==0 else 0
        elif key==ord('r'): #重新开始
            return 1
    return 0

while runWorld(100,100,c=20,wait=0):
    pass

input("Press Enter to Continue")

cv2.destroyAllWindows()